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Abstract 

We redraw, using state-of-the-art methods for free-energy calculations, the phase diagrams of 
two reference models for the liquid state: the Gaussian and inverse-power-law repulsive potentials. 
Notwithstanding the different behavior of the two potentials for vanishing interparticle distances, 
their thermodynamic properties are similar in a range of densities and temperatures, being ruled by 
the competition between the body-centered-cubic (BCC) and face-centered-cubic (FCC) crystalline 
structures and the fluid phase. We confirm the existence of a reentrant BCC phase in the phase 
diagram of the Gaussian-core model, just above the triple point . We also trace the BCC-FCC 
coexistence line of the inverse-power-law model as a function of the power exponent n and relate the 
common features in the phase diagrams of such systems to the softness degree of the interaction. 
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I. INTRODUCTION 



A recurrent issue in statistical physics is the prediction of the phase behavior of a macro- 
scopic system from the knowledge of the force law which rules the interaction between the 
constituent particles. Over the last decades, significant progress has been made in fulfilling 
this task, at least for simple model systems of classical point particles interacting through a 
spherically-symmetric pair potential v(r). In particular, continuous advancements in both 
computer simulation methods and density-functional theory have revealed a few general fea- 
tures of the freezing transition, a phase change that is prominently driven by the repulsive 
component of the interatomic potential. For assigned thermodynamic conditions, the soft- 
ness of the repulsion substantially determines the nature of the incoming solid phase. In fact, 
a rigid potential favors a more compact phase, such as the face-centered-cubic (FCC) or the 
hexagonal-close-packed (HCP) crystalline structure, as opposed to, say, the body-centered- 
cubic (BCC) one. This is largely true even for more complex substances such as those that 
are generically classified as soft materials, like colloidal suspensions for which the effective 
pair interactions can be accurately tailored through a fine experimental control of suspended 
particle and solvent properties ^ Q| . A scholastic example is offered by the entropy-driven 
repulsion between (the centers of mass of) self-avoiding polymer coils dispersed in a good 
solvent 

Baa 

. In this specific case, the repulsion happens to be finite even if "particles" 
fully overlap and can be modelled with a Gaussian-shaped potential 



In Eq. (1.1) the length a is of the order of the gyration radius of the polymer while the 
(positive) energy e is of the order of k^T, where T is the temperature and k-g is the Boltzmann 
constant. However, it is also useful to study this potential without assuming an explicit 
dependence of e and a upon the temperature. The model is then referred to as the Gaussian- 
core model (GCM) [5]. Other examples of softly repulsive fluids are provided by dilute 
solutions of colloidal particles, sterically stabilized against coagulation through a polymer 
coating. In this case, the strength and range of the effective repulsion between a pair of 
particles are determined by the thickness of the layer of grafted polymer chains Q and can 
be modelled with an inverse-power-law (IPL) potential 




(1.1) 




(1.2) 
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upon choosing a suitable power exponent n. In the past, this potential has been also used 
for a simplified description of the thermal behavior of various simple metals under extreme 
thermodynamic conditions 0, 

With the advent of numerically accurate techniques for the calculation of the Helmholtz 
free energy of a solid, computer simulation has become a valuable tool for the determination 
of fluid-solid and solid-solid coexistence lines. Most notably, the Frenkel-Ladd (FL) method 
can be used in conjunction with conventional thermodynamic integration techniques to cal- 
culate the "exact" phase diagram of the system, the only source of error being ultimately 
associated with the quality of the statistical sampling of the (finite) system 
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fledged computer studies have been carried out for the repulsive Yukawa potential 



. Full- 
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141, and for the 



17j. In this 



typically used to model solutions of charged colloids with counterions 
logarithmic-plus- Yukawa potential modelling star-polymer interactions 
paper, we present accurate results for the phase diagrams of the Gaussian-core model and of 
the inverse-power-law potential, with the aim of providing a comparable benchmark for such 
two potentials. A preliminary report on the present findings for the GCM has been already 
given in [18]. The main results of this study are: i) the discovery of a reentrant BCC phase in 
the phase diagram of the GCM; ii) an accurate determination of the BCC-FCC coexistence 
point of the IPL potential as a function of the softness parameter 1/n; iii) the formulation of 
an argument that accounts for the similarity of the low-density/low-temperature behavior 
exhibited by the two models. 



II. MODELS AND METHOD 



A. The Gaussian-core model 



As is well known, the freezing transition of a many-particle system is ultimately pro- 
moted by the Pauli repulsion between inner-shell electrons, causing the effective interatomic 
potential to diverge at short distances. It is probably less known that the existence of a 
thermodynamically stable solid does not require a singular repulsion for vanishing inter- 



atomic separations. In fact, even a finite square- 
support a stable solid phase at all temperatures 



carrier potential (penetrable spheres) can 
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20j . In this respect, the Gaussian po- 



tential is a more realistic form of bounded repulsion. The GCM was originally introduced 
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by Stillinger [sj. Two distinctive features of the phase diagram of the GCM, that are absent 
in other bounded-interaction models, are the existence of a maximum melting temperature, 
T max , and, correspondingly, of a reentrant melting transition into a denser fluid phase for 



T < T 

1 ~- * TT 
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231 ] . In his original paper on the GCM, Stillinger noted that, in the limit 



of vanishing temperature and density, the Gaussian particles behave as hard spheres with a 
larger and larger diameter. In the same limit, the fluid freezes into a FCC structure whose 
number density p vanishes with the temperature T (from now on, temperature, pressure, and 
density will be given in reduced units, i.e., T* = k B T/e, P* = Pcr 3 /e, and p* = per 3 ). The 
calculation of the total energy of different cubic structures shows that the FCC structure is 
favored, at zero temperature, for reduced densities lower than 7r -3 / 2 ~ 0.1796. Beyond this 
threshold, a BCC solid takes over. However, upon compression, any regular arr ang ement of 



particles is eventually doomed to collapse for any T > (reentrant melting) |2lJ, |24j . An 
extended study of the GCM phase diagram was carried out by Lang and coworkers who 
employed an integral-equation theory for the fluid phase and a harmonic approximation 
for the crystalline phases j^ . In this combined approximation, these authors estimated a 
maximum melting temperature T^ ax = 0.0102 at which the coexisting fluid and BCC phases 
have a reduced density p^ax = 0.2292. They also found that the FCC-BCC coexistence lines 
run linearly from the points Pfcc = 0.17941 and Pbcc = 0.17977 at T* = to the points 
Pfcc = 0.16631 and p^cc = 0.16667 at the triple-point temperature T t * = 0.008753. 



B. The inverse-power-law fluid 



The IPL potential has been the subject of many studies [25J, |26|, |27(. This potential 
provides a continuous path from hard spheres (n — > oo) to the one-component plasma 
(n = 1). Actually, a neutralizing background is needed for n < 3 to balance the non- 
integrable long-ranged repulsion J3, |3| . The thermodynamic properties of the model can 
be expressed in terms of just one quantity, 7 = p*T*~ 3 / n , as a result of the scaling properties 
of the excess Helmholtz free energy. Hence, for any given n, it is sufficient to trace the thermal 
behavior of the IPL fluid as a function of p* along the isotherm T* — 1. For large values 
of n, the BCC phase is found to be unstable with respect to shear modes. As a result, the 
fluid freezes into a FCC solid. As n decreases, the potential becomes increasingly softer and 
longer ranged until, for n w 8, the BCC phase becomes mechanically stable. For smaller 
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values of n and for densities close to freezing, the entropy of the BCC solid turns out to be 
greater than that of the FCC phase since the former structure can actually sustain a larger 
number of soft shear modes. However, the FCC phase is favored at high densities and low 
temperatures because of its smaller energy. Hence, upon reducing 7, a FCC-BCC transition 
may actually anticipate the melting of the solid for intermediate values of n. Agrawal and 
Kofke located the triple point in the range 0.11 < 1/n < 0.17, further claiming that a more 
realistic estimate of the lower bound, based on the mechanic stability range of the BCC 
phase, was 0.1305 (2?}]. However, they did not provide a direct and systematic estimate of 
the location of the FCC-BCC transition as a function of the power exponent n. 

C. Details of the Monte Carlo simulation 

We performed Monte Carlo (MC) simulations of the GCM and of the IPL model in the 
canonical ensemble, implementing the standard Metropolis algorithm with periodic bound- 
ary conditions and with the nearest-image convention. The numbers of particles N was 
chosen so as to fit a cubic box of volume V with an integer number of cells: iV = 4m 3 for 
the FCC solid and N = 2m 3 for the BCC solid, m being the number of cells along any 
spatial direction. The nearest-neighbor distance a is (V2/2)(L/m) for a FCC lattice and 
(y/3/2)(L/m) for a BCC lattice, where L = (iV/p) 1 / 3 . Our typical samples consisted of 1372 
particles in the fluid phase and in the solid FCC phase, and of 1458 particles in the BCC 
phase. We emphasize that, unless one resorts to some extrapolation to infinite size (a com- 
putationally demanding task), one should compare the statistical properties of samples with 
similar sizes, in particular if one wants to investigate the relative thermodynamic stability 
of different phases at given temperature and pressure. Occasionally, we considered smaller 
as well as larger sizes to check whether our results could be possibly affected by a significant 
size dependence (and they were actually not). We also carried out a number of simulations 
of a HCP solid hosting 10 x 12 x 12 = 1440 particles. 

In order to locate the freezing point at a given temperature, we generated distinct chains 
of simulations, starting from the very dilute fluid on one side and from the dense solid on 
the other side. As a rule, the last MC configuration produced for a given value of p was 
used, after a suitable rescaling of particle coordinates, as the starting point for the next run 
at a slightly higher or lower density. We monitored the difference in energy and pressure 
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between the various phases until we observed a sudden change of both quantities, so as 
to avoid averaging over heterogeneous thermodynamic states. Usually, the solid can be 
overheated for a few more steps beyond the melting point. The undercooling of the fluid is 
notoriously much easier. For any p and T, the equilibration of each sample typically took 
a few thousand MC sweeps, a sweep consisting of one attempt to change sequentially the 
position of all the particles. The maximum random displacement of a particle in a trial MC 
move was adjusted once a sweep during the run so as to keep the acceptance ratio as close 
to 50% as possible, only small excursions around this value being allowed. Thermodynamic 
averages were computed over trajectories 3 x 10 4 sweeps long. The pressure was computed 
through the virial formula 



where is the distance between particles % and j and (• • • ) denotes the canonical aver- 
age. For the IPL fluid, the excess pressure is proportional to the excess energy u cx since 
—rv' n (r) = nv n (r). To avoid counting pair interactions twice, the potential was truncated at 
a cutoff distance r c , which was taken to be just slightly smaller than L/2. Correspondingly, 
long-range corrections to the energy and pressure were calculated upon assuming the radial 
distribution function (RDF) g(r) to be unity for r > r c . The histogram of g{r) was con- 
structed with a spatial resolution Ar = cr/50 and was updated every 10 MC sweeps. The 
RDF was never found to differ significantly from unity - in any phase - at the maximum 
distance L/2, at least for the largest sample sizes. 

The numerical estimate of statistical errors is a critical issue whenever different candidate 
solid structures so closely compete for thermodynamic stability. To this aim, we divided the 
MC trajectory into ten blocks and estimated the error bars to be twice as large as the stan- 
dard deviation of the block averages, under the implicit assumption that the decorrelation 
time of any relevant variable was less than the size of a block. Tipically, the relative errors 
affecting the energy and pressure of the fluid were found to be a few hundredths percent 
at the most (for a solid, they were even smaller). We computed the difference between the 
excess free energies of any two equilibrium states (say, 1 and 2) belonging to the same phase 
with the standard thermodynamic integration technique, i.e., via the combined use of the 




(2.1) 
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formulae 

,/ex(P, T 2 ) /ex(p, Ti) /" T2 M cx (p, T) 



and 



fp2 1 

l;U f >2-T) = iM, H -T)+ I dp- 



pi 



PP(p,T) 
P 



(2.3) 



In the above equations, / ex is the excess Helmholtz free energy per particle and f3 = (k B T)~ l . 
The integrals in Eqs. (2.2) and (2.3) were performed numerically with Simpson's rule over a 
cubic-spline approximant of the simulation data. The excess values of the chemical potential 
and of the entropy per particle follow from 

/3/iex = Pfex + — - 1 and ^ = /3(w ex - / ex ) . (2.4) 
P k B 

Equations (2.2) and (2.3) are useless if an independent estimate of the free energy is not 
available for some selected reference state of each phase. The free energy of a solid can be 
obtained with the FL method, that is reviewed in the next paragraph. As for the fluid, a 
reference state can be any equilibrium state characterized by a very small density (a nearly 
ideal gas), since then the excess chemical potential be accurately estimated through 



Widom's particle-insertion method 
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/i ex = -k B T In (exp (-(3E ixm )) , (2.5) 

where E ins includes all the interactions between a randomly inserted test particle and the 
particles of the fluid. Typically, the average in Eq. (2.5) was computed over equilibrium 
runs 6 x 10 4 sweeps long, an insertion being attempted at the completion of each sweep. 

The chemical potential is the quantity that decides which phase is stable for given values 
of T and P. Hence, it is necessary to inspect the error affecting this quantity with care. 
Let us consider, for instance, the chemical potential of a solid whose density is p 2 in a state 
that is connected, along an isothermal path, to another (reference) state with density p\. 
Let Sfi be the numerical precision of the FL calculation of the free energy of the reference 
solid. For each NVT simulation carried out along the isothermal path, one can estimate the 
error affecting the compressibility factor (3P/ p, whence its mean error 5z along the path. 
Hence, according to Eqs. (2.3) and (2.4), the error affecting p is approximately given by 
8f 1 + k B T5z(l + Hp 2 / Pl )). 
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D. The Frenkel-Ladd calculation of the solid free energies 



The method proposed by Frenkel and Ladd to calculate the free energy of a solid relies 
on a different kind of thermodynamic integration Q|. The idea is to connect continuously 
the state of a solid (say, system 1 with potential energy Vi), to the state of another solid 
(say, system with potential energy Vq) whose free energy is known. Upon assuming the 
reference solid to be an Einstein solid with spring constant c, a hybrid system is sampled 
with potential V\ = Vq + \{Vi — Vq) (0 < A < 1), for given p and T. The Helmholtz free 
energy of system 1 is then obtained through Kirkwood's formula: 

F 1 -F Q = [ dA (V, - V Q ) X . (2.6) 
Jo 

The free energy of the Einstein solid and the mean square separation of an Einstein particle 
from its reference lattice site are: 

AJ$ = |;, (2-8) 

where A is the thermal wavelength. A nontrivial problem with this choice of Vq, already 
pointed out in the original paper 8|, is the following: The Einstein particles can only perform 
limited excursions about their reference lattice positions, which implies a finite value of AR 2 ; 
however, there is no means to constrain particles interacting with the potential V\ to move 
within the neighborhood of their initial positions, even in the solid phase. In other words, at 
variance with Vo, V\ is a translationally-invariant potential, which implies that the integrand 
in Eq. (2.6) diverges for A = 1. This problem can be overcome by performing MC simulations 
of V\ under the constraint of a fixed center of mass. Only afterwards, the proper corrections 
are taken into account We just quote the final result: 

P /",, l„ „<o) ,A CM 



where the superscript CM denotes a constrained average and is the total potential 
energy V\ for particles located in their respective lattice positions. A linear scaling relation 
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was conjuctered, on rather general grounds, for the asymptotic behavior of the free energy 
as a function of the sample size: (3f ex {N) + In N/N = (3f ex {oo) + C^iV" 1 ) Q. 

We finally add some further considerations about the numerical implementation of the 
FL method. In principle, one might choose the value of c in an arbitrary way. However, if 
the quantity AR 2 for the target solid is close to 3/ (/3c), the variations of the integrand in 
Eq. (2.9) and, correspondingly, the error made in the numerical quadrature get very small. 
In practice, A was taken to increase from to 0.9 with steps of 0.05, and with a smaller step 
(0.01) in the range from 0.9 to 1. For each value of A, as many as 2 x 10 4 MC sweeps were 
produced at equilibrium. The MC moves were such that the position of the center of mass 
of the solid was held fixed; this way, particles can undergo only limited excursions along the 
run. For this reason, if a particle happens to move out of the simulation box, there is no 
need to put it back. 



III. RESULTS 



A. Gaussian-core model 



Figure 1 presents the phase diagram of the GCM, calculated with the numerical tech- 
niques discussed above, in the temperature-pressure plane. The general features of this 
phase diagram have already been illustrated in a number of papers |3, [3, [yj ^ • A novel 
feature that has been just reported is the emergence, at low densities and in a narrow range 
of temperatures (0.0031 < T* < 0.0037), of a BCC phase intermediate between the fluid 
and the FCC phase Correspondingly, the thermodynamic coordinates of the triple 

point are: T t * = 0.0031 and P* T = 0.0070. Calculations carried out with different number of 
particles confirm that the reentrant BCC phase at low densities is not a finite-size effect (see 
the inset in Fig. 1). Tracing the coexistence lines in a phase diagram requires the calculation 
of the chemical potential, fi(T,P), for all candidate phases. In the case of the GCM, the 
reentrance of the fluid phase at high densities for T < T max , calls for an accurate calcula- 
tion of this quantity in the dense-fluid regime. However, in such thermodynamic conditions 
Widom's method is no longer reliable. We thus resorted to a thermodynamic integration 
along a composite path in the fluid-phase region: We first moved from low to high densities 
(p* = 0.74) at fixed temperature (T* = 0.015), and followed therefrom an isochoric path 
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down to low temperatures (T* = 0.002). Differences between the chemical potentials of fluid 
and solid phases were plotted as a function of the pressure in Figs. 2 and 3. The sequence 
of phases showing up with increasing pressures at T* = 0.002 is fluid, FCC, BCC, fluid 
again. It is also apparent from Fig. 2 that the chemical-potential gap between the FCC and 
BCC phases depends on the pressure in a distinctly nonmonotonic way: At low pressures, 
the BCC phase already competes with the FCC phase and is about to become stable. The 
changeover from the FCC to the BCC phase occurs at higher temperatures. In fact, the 
shallow valley in the quantity /ifcc — ^bcc gradually slides upwards with the temperature 
until, for T* > 0.0030, the BCC phase slips in between the fluid and the FCC phases. As 
also seen from Fig. 3, the HCP phase never comes into play. 

The excess Helmholtz free energies calculated for some solid phases of the GCM are given 
in Table 1, together with the values of the reduced elastic constant c* = ccr 2 /e. Though one 
single FL calculation would be sufficient to estimate the free energy of any given phase in 
any other thermodynamic state, we found it more practical to repeat the calculation rather 
than generating a dense mesh of thermodynamic paths. As illustrated in Fig. 4, the free 
energy scales asymptotically with N as conjectured by Poison and coworkers [9j. 

The thermodynamic properties of the GCM at fluid-solid and solid-solid coexistence 
are given in Table 2. Upon increasing the pressure for temperatures just above the triple 
point, the entropy per particle drops by about 0.8 k-Q and 0.1 k-Q across the fluid-to-BCC and 
BCC-to-FCC phase transitions, respectively. However, the next jump at the FCC-to-BCC 
transition is positive in that the entropy of the coexisting BCC phase, notwithstanding its 
slightly higher density (see Table 2), again exceeds that of the FCC solid by about the 
same amount (0.1 k-o). As then implied by the Clausius-Clapeyron equation, the slope of 
the higher-pressure branch of the FCC-BCC coexistence line is negative in the (T, P) plane 
as is also the slope of the upper branch of the BCC- fluid coexistence line. In fact, at the 
reentrant melting point, the coexisting fluid has both a larger density and a larger entropy 
than the BCC phase, the entropy jump being again of the order of 0.8 k-Q per particle. 

We actually found that for low densities (p* < 0.2) the entropy of the BCC solid system- 
atically overcomes that of the FCC solid, even below the triple-point temperature. The gap 
decreases slowly with increasing temperatures and is likely due to the larger number of soft 
shear modes hosted by the BCC structure. Actually, even at low temperatures, BCC-ordered 
grains are expected to form in the freezing fluid, a circumstance which may slow down the 
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crystallization process. On the energy side, the scales tip instead in favor of the FCC solid. 
In fact, at low densities (p* < 0.17), the FCC phase has a lower energy than the BCC phase 
for T* < 0.008. Not surprisingly, this temperature is very close to the reduced triple-point 
temperature predicted by Lang and coworkers 2^. These authors actually used the Gibbs- 
Bogoliubov inequality to optimize a harmonic model of both solid phases. Considering the 
fact that the FCC-BCC transition occurs at rather low densities, the above approximation 
does likely underestimate the (relative) entropies of the competing crystalline phases. As 
a result, the FCC phase turns out to be stable in a much more extended region than that 
predicted by free-energy calculations carried out with the FL method (see also Fig. 2 of Q|). 

As shown in Fig. 1, we predict that the BCC phase is stable for temperatures lower than 
^max — 0.00874. The reduced density and pressure of the solid coexisting with the fluid at 
this temperature are p^ ax ~ 0.239 and -P^ax — 0.128, respectively. At this bending point 
the BCC solid melts, upon heating, with no volume change. However, the transition is still 
first-order since the entropy per particle of the fluid is 0.79 /cb larger than in the solid. As for 
the FCC phase, our current estimate for the uppermost FCC-BCC coexistence temperature 
is 0.0038. The corresponding estimates for the reduced density and pressure attained by the 
system in this extremal state of the FCC stability basin are 0.123 and 0.0174, respectively. 
Comments analogous to those made above on the nature of the BCC-fluid phase transition 
at the maximum temperature of the coexistence line can be straightly extended to the 
FCC-BCC transition as well. 

Figure 5 shows the RDFs of the fluid and solid phases for T* = 0.004, at both ordinary and 
reentrant melting densities. We note that neighboring particles lie, on average, at relative 
distances larger than ct/a/2, which is where the Gaussian potential inflects. 



B. Inverse-power-law model 

We list the excess Helmholtz free energies for some crystalline states of the IPL model in 
Table 3. Figure 6 shows the chemical potential plotted as a function of the pressure for an 
intermediate value of the softness parameter, 1/n = 0.16. Upon increasing 7, the IPL fluid 
freezes into a BCC structure. A further compression triggers a solid-solid transition into a 
FCC structure. Correspondingly, the equation of state (Fig. 7) shows two tie lines joining 
the three distinct branches at the fluid-BCC and BCC-FCC transition points. We found 
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that the HCP solid is systematically less stable than the FCC solid, although the chemical- 
potential gap between the two phases is rather small at all pressures. As n grows, the range 
of BCC stability gradually reduces until, for 1/n < 0.14, the FCC structure only is left over. 
The maximum value of the power exponent for which the BCC solid is thermodynamically 
stable is close to 7 (0.14 < 1/n < 0.15). The phase diagram of the IPL fluid is shown in 
Fig. 8 while the phase-transition densities are listed in Table 4. In the triple-point region, 
the density jumps are approximately 0.030 across the ffuid-BCC transition and 0.004 across 
the BCC-FCC transition; correspondingly, the entropy decreases at each transition by about 
0.7 fee and 0.1 /cb per particle, respectively. All along the freezing line, the entropy of the 
BCC solid overcomes that of the FCC solid by approximately 0.1 k-Q per particle. We verified 
that the estimate of the power exponent n tT for which the IPL model exhibits triple-phase 
coexistence is not affected by appreciable finite-size effects because of the large samples 
used, the only significant source of error being the limited statistics of the simulation. The 
maximum error expected on /3A/i, near to freezing, is approximately equal to 1.5 x 10~ 2 , 
as for the GCM |l£|. Such an error would, in turn, affect the estimates of the reduced 
phase-transition densities by an amount not greater than 3 x 10~ 2 . We conclude that the 
error on n tr should not be larger than 0.5. 

Figure 8 presents the phase diagram of the IPL model, spanned by the softness parame- 
ter 1/n and by the scaled quantity 7. We also plotted fluid-solid coexistence data obtained 



by other authors using different approaches |2i 



3 n 



2fiL I27L l33| . The data reported in |25j refer 



to a single-occupancy solid in which each particle is confined within a spherical cell whose 
diameter is set equal to the FCC nearest-neighbor lattice spacing. Dubin and Dewitt esti- 
mated the fluid-solid and solid-solid coexistence of the IPL model by applying the empirical 
Lindemann criterion, including some anharmonic contributions to the internal energy of the 
solids (34}]. Their fluid-FCC coexistence line almost coincides with the present fluid-solid 
coexistence line; instead, the location of the BCC threshold does not agree with the present 
data. 



IV. DISCUSSION 

The low-density/low-temperature topology of the phase diagram of the GCM, with a 
triple point separating a region where the fluid freezes into a FCC structure from another 
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region where such two phases are bridged by an intermediate BCC phase, is common to 
other model systems with softly repulsive interactions. Another example, besides the GCM 
and the IPL model, is the Yukawa potential, vy(r) = Bexp(—r/£)/r (with B > 0), which 
has been used to model the thermal behavior of dilute solutions of charged colloids with 



counterions 
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131 ] . The phase diagram of this latter model can be conveniently spanned 



by two scaled quantities, i.e., p£ 3 and Bp 1 ! 3 / '{k^T) [l2j. In an attempt to relate the phase 
behaviors of the Gaussian and Yukawa potentials to that of the IPL model, we required in 
the logarithmic derivatives of the three potentials to match, at least for separations close to 
the average nearest-neighbor distance, r = p~ 1//3 . This condition yields two "optimal" values 
for the IPL power exponent: uq = 2(pcr 3 ) _2 / 3 and ny = l+(p£ 3 )~ 1//3 , respectively. Increasing 
the density of the Gaussian and Yukawa models along the freezing line is tantamount, 
through the above relations, to an increase of the softness parameter 1/n of the "equivalent" 
IPL model. In both cases the system is driven towards the BCC phase. The values of uq 
and calculated at the triple point of the Gaussian and Yukawa potentials are 9.6 and 
12.0, respectively, not too far from the IPL value 7). This rough argument shows that it 
is possible to gauge the softness of such potentials in such a way that their phase behavior is 
referred to a common frame where the BCC phase is ultimately promoted by a sufficiently 
soft interaction. 

A similar interplay between BCC and FCC solid phases is also present in the phase dia- 
gram of the modified Buckingham potential that is used to model the behavior of rare gases 
at high temperature and pressure js^ : under extreme thermodynamic conditions close to 
freezing, notwithstanding the presence of an attractive tail, particles feel almost exclusively 
the repulsive shoulder of the potential, whose softness increases with the temperature. 



V. CONCLUSIONS 



In this paper we revisited the phase diagram of two classical reference models for the liquid 
state, both governed by softly repulsive interactions with Gaussian and inverse-power-law 
shapes, respectively. The aim of this effort was to point out some aspects of the phase 
diagrams which had been previously overlooked or simply ignored, in a common theoret- 
ical framework for repulsive pair potentials that decrease polynomially or faster with the 
interparticle distance. To this end, we carried out extensive Monte Carlo simulations and 



12 



calculated the thermodynamic properties of both fluid and solid phases through the com- 
bined use of thermodynamic integration techniques and free-energy calculations. We further 
discussed the similarities and semi-quantitative correspondences in the phase diagrams of 
such models in relation to the softness of the potential and with specific reference to the 
interplay of fluid, BCC and FCC crystalline phases. 
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TABLE CAPTIONS 



Table 1 : Excess Helmholtz free energy per particle / ex , in units of JcbT, for some crystalline 
states of the GCM: FCC (N = 1372), HCP (N = 1440), and BCC (N = 1458). For 
T* = 0.003 and 0.006, the tabulated values refer to systems with 864 and 1024 particles, 
respectively. The value of the reduced elastic constant c* = ccr 2 /e, playing a role in 
the Frenkel-Ladd calculation, is also displayed in square brackets: for such values of c, 
the mean square displacement of the Einstein solid approximately matches the mean 
square deviation of a GCM particle from its reference lattice site. 

Table 2 : Thermodynamic properties of the GCM at fluid-solid and solid-solid coexistence. 
With the only exception of the states identified by T* = 0.003 and 0.006, the simulation 
data refer to samples with N = 1372 (fluid and FCC) and N = 1458 particles (BCC). 
For T* = 0.003 and 0.006, the tabulated values are relative to samples with 864 and 
1024 particles, respectively. We show, from left to right, the values relative to the 
following first-order phase transitions: freezing of the fluid into either a FCC or BCC 
solid; BCC-to-FCC transition; FCC-to-BCC transition, and high-density melting of 
the BCC solid. We indicated with p f and p m the freezing and melting densities, 
respectively. All data were arbitrarily rounded off; note, however, that the error 
originating from the limited precision of both the Monte Carlo data and the free-energy 
values may be actually larger than that apparently implied by the chosen truncation. 

Table 3 : Excess Helmholtz free energy per particle / ex , in units of A; B T, for some crystalline 
states of the IPL model: FCC (N = 1372), HCP (N = 1440), and BCC (N = 1458). 
The value of the reduced elastic constant c* = ca 2 /e is also displayed in square brackets 
(see the caption of Table 1). 

Table 4 : Thermodynamic properties of the IPL model at fluid-solid and solid-solid coex- 
istence for some values of the softness parameter 1/n. The tabulated data refer to 
samples with N = 1372 (fluid, FCC) and N = 1458 particles (BCC). All data were 
arbitrarily rounded off at the third decimal digit (see the caption of Table 2). 
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FIGURE CAPTIONS 



Fig. 1 : Phase diagram of the GCM in the (T, P) plane: The markers (open circles) identify 
the states were we carried out the numerical calculations with N = 1372 particles for 
the fluid and FCC phases and with N = 1458 for the BCC phase. The triple-point 
region is magnified in the inset: The solid circles mark the phase thresholds in larger 
samples with 2048 particles for the fluid and FCC phases and with 2000 particles for 
the BCC phase, at T* = 0.0037. The FCC-BCC coexistence pressure at T* = (open 
square) was taken from |5]. The lines drawn through the data points are a guide for 
the eye. 

Fig. 2 : Chemical-potential differences between pairs of GCM phases plotted as a function 
of the pressure for T* = 0.002: /9A/Xfl U id,Fcc (dashed line); /9A/ifl ui d,Bcc (dotted line); 
/3A/ipcc,BCC (continuous line). Upon increasing either P or p, the GCM exploits a 
fluid-FCC-BCC-fluid sequence of phases. A zoom in the low-pressure region reveals 
the non-monotonic behavior of /3Ap F cc,BCC; a feature that is ultimately responsible 
for the reentrance of the BCC phase at higher temperatures (see also Fig. 3). The 
lines are spline interpolants of the data. 

Fig. 3 : Chemical potential differences between pairs of GCM phases plotted as a function 
of the pressure for T* = 0.0037: /3Apfl ui d,Fcc (dashed line), /3Apfl ui d,Bcc (dotted line), 
/3Apfcc,hcp (dash-dotted line), and /3A/zfcc,bcc (continuous line). The lines for 
/3A/ifl U id j fcc an d /3AyU fluid B cc are hardly resolved on the scale of the figure. Upon 
increasing the pressure, the GCM exploits a fluid-BCC-FCC-BCC-fluid sequence of 
phases. The lines are spline interpolants of the data. 

Fig. 4 : Frenkel-Ladd calculation of the free energy of the GCM solid. Top: The integrand 
of Eq. (2.9) plotted as a function of A for a BCC solid with 1458 particles (p* = 0.3, 
T* = 0.004); the continuous line is a spline interpolant of the data, the error bars being 
much smaller than the size of the symbols. Bottom: The quantity /3f ex (N) + In N/N 
plotted as a function of N" 1 for a BCC solid in the same thermodynamic state, with 
iV = 432, 686, 1024, 1458. We checked the asymptotically linear scaling of the above 
quantity also for other states among those listed in Table 1. 
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Fig. 5 : Radial distribution functions of fluid and solid phases of the GCM for T* = 0.004, 
p* = 0.11 (top) and p* = 0.52 (bottom): fluid (continuous line), BCC (dotted line), 
FCC (dashed line). 

Fig. 6 : Chemical-potential differences between pairs of phases of the IPL fluid for l/n — 
0.16, plotted as a function of the pressure for T* = 1: /3A/z flu id,Fcc (dashed line); 
/^A/iflui^BCC (dotted line); /?A/ifcc,hcp (dash-dotted line); /3A/ifcc,bcc (continuous 
line). For increasing pressures, the IPL fluid exploits a ffuid-BCC-FCC sequence of 
phases. The lines are spline interpolants of the data. 

Fig. 7 : Equation of state of the IPL model for T* = 1 and l/n = 0.16. The fluid branch 
(O) continuous line) and the FCC branch (□, dashed line) were plotted for iV = 1372 
particles, while the BCC branch (A, dotted line) was obtained with a sample of 1458 
particles. The coexistence densities, signalled by two pairs of thin vertical lines, are 
located where the tie lines cross different branches of the equation of state. 

Fig. 8 : Phase diagram of the IPL potential: fluid (N = 1372, open circles); FCC 
(N = 1372, open squares); BCC (N = 1458, open triangles). The continuous lines 
drawn through the data are a guide for the eye. The dashed lines are the phase 



thresholds reported in 



331 ] . In the inset, we compare the present data (continuous 



lines) with the phase-transition loci reported in ^5l (n = 6: fluid, solid circle; FCC, 

n n 

solid square); in 26] {n = 6: fluid, star; BCC, tripod; FCC, cross); and in [2j| (fluid, 
cross; solid, tripod; note that this latter datum refers to FCC for l/n < 0.16 and to 
BCC otherwise). 
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